EMODnetWFS Case Study: Accessing and mapping EMODnet data
Source:vignettes/emodnetwfs.Rmd
emodnetwfs.RmdIntroduction
The package was designed to make EMODnet vector data layers easily accessible in R. The package allows users to query information on and download data from all available EMODnet Web Feature Service (WFS) endpoints directly into their R working environment. Data are managed as sf objects which are currently the state-of-the-art in handling of vector spatial data in R. The package also allows user to specify the coordinate reference system of imported data.
Data Product
Installation
You can install the development version of EMODnetWFS from GitHub with:
remotes::install_github("EMODnet/EMODnetWFS")Explore the EMODnet WFS services with R
For this tutorial we will make use of the sf, dplyr and mapview packages. The simple features sf package is a well known standard for dealing with geospatial vector data. The package dplyr is a strong library for data manipulation. This package also loads magickr’s pipe operator %>%, which allows to write pipelines in R. To visualize geometries, mapview will create quick interactive maps.
Run this line to install these packages:
install.packages(c("sf", "dplyr", "mapview"))With the EMODnetWFS package, we can explore and combine the data served by the EMODnet lots through OGC Web Feature Services or WFS.
Imagine we are interested on seabed substrates. The first stop is to choose what EMODnet lot can provide with these data. For that, we can check the services available with the emodnet_wfs() function.
library(EMODnetWFS)
library(mapview)
library(dplyr)
library(sf)
emodnet_wfs()
#> # A tibble: 17 × 2
#> service_name service_url
#> <chr> <chr>
#> 1 bathymetry https://ows.…
#> 2 biology https://geo.…
#> 3 biology_occurrence_data https://geo.…
#> 4 chemistry_cdi_data_discovery_and_access_service https://geo-…
#> 5 chemistry_cdi_distribution_observations_per_category_and_region https://geo-…
#> 6 chemistry_contaminants https://nodc…
#> 7 chemistry_marine_litter https://www.…
#> 8 geology_coastal_behavior https://driv…
#> 9 geology_events_and_probabilities https://driv…
#> 10 geology_marine_minerals https://driv…
#> 11 geology_sea_floor_bedrock https://driv…
#> 12 geology_seabed_substrate_maps https://driv…
#> 13 geology_submerged_landscapes https://driv…
#> 14 human_activities https://ows.…
#> 15 physics https://geos…
#> 16 seabed_habitats_general_datasets_and_products https://ows.…
#> 17 seabed_habitats_individual_habitat_map_and_model_datasets https://ows.…The column service_name shows services available, while service_url has the corresponding base url to perform a WFS request. The Seabed portal should have the data we are looking for. A WFS client can be created by passing the corresponding service_name to the function emodnet_init_wfs_client(). The layers available to this WFS client are consulted with emodnet_get_wfs_info().
seabed_wfs_client <- emodnet_init_wfs_client(service = "seabed_habitats_general_datasets_and_products")
#> Loading ISO 19139 XML schemas...
#> Loading ISO 19115 codelists...
#> Loading IANA mime types...
#> No encoding supplied: defaulting to UTF-8.
#> ✔ WFS client created succesfully
#> ℹ Service: 'https://ows.emodnet-seabedhabitats.eu/geoserver/emodnet_open/wfs'
#> ℹ Version: '2.0.0'
emodnet_get_wfs_info(wfs = seabed_wfs_client)
#> # A tibble: 33 × 9
#> data_source service_name service_url layer_namespace layer_name title
#> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open art17_hab… 2013…
#> 2 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open art17_hab… 2013…
#> 3 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open art17_hab… 2013…
#> 4 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open art17_hab… 2013…
#> 5 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open art17_hab… 2013…
#> 6 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open art17_hab… 2013…
#> 7 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open art17_hab… 2013…
#> 8 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open art17_hab… 2013…
#> 9 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open biogenic_… Biog…
#> 10 emodnet_wfs seabed_habitats_gen… https://ow… emodnet_open habitat_p… Coll…
#> # … with 23 more rows, and 3 more variables: abstract <chr>, class <chr>,
#> # format <chr>Each layer is explained in the abstract column. We can see several layers with the information provided by the EU member states for the Habitats Directive 92/43/EEC reporting. We will select the layers about coastal lagoons, mudflats and sandbanks with their respective layer_name.
habitats_directive_layer_names <- c("art17_hab_1110", "art17_hab_1140", "art17_hab_1150")
emodnet_get_layer_info(wfs = seabed_wfs_client, layers = habitats_directive_layer_names)
#> # A tibble: 3 × 9
#> data_source service_name service_url layer_namespace layer_name title abstract
#> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 emodnet_wfs https://ows… seabed_hab… emodnet_open art17_hab… 2013… "Gridde…
#> 2 emodnet_wfs https://ows… seabed_hab… emodnet_open art17_hab… 2013… "Gridde…
#> 3 emodnet_wfs https://ows… seabed_hab… emodnet_open art17_hab… 2013… "Gridde…
#> # … with 2 more variables: class <chr>, format <chr>We are now ready to read the layers into R with emodnet_get_layers(). EMODnetWFS reads the geometries as simple features (See sf package) transformed to 4326 by default. Specifying another map projection is possible by passing a EPGS code or projection string with emodnet_get_layers(crs = "your projection"). The argument reduce_layers = TRUE stack all the layers in one single tibble. Default is FALSE and returns a list of sf objects, one per layer.
habitats_directive_layers <- emodnet_get_layers(wfs = seabed_wfs_client,
layers = habitats_directive_layer_names,
reduce_layers = TRUE)
class(habitats_directive_layers)
#> [1] "sf" "data.frame"
glimpse(habitats_directive_layers)
#> Rows: 221
#> Columns: 9
#> $ gml_id <chr> "art17_hab_1110.13", "art17_hab_1110.22", "art17_h…
#> $ habitat_code <chr> "1110", "1110", "1110", "1110", "1110", "1110", "1…
#> $ ms <chr> "DK", "ES", "ES", "PT", "PT", "PL", "DK", "FR", "U…
#> $ region <chr> "ATL", "MAC", "MMAC", "MMAC", "MATL", "MBAL", "MBA…
#> $ cs_ms <chr> "U2+", "U1+", "U1+", "XX", "U1-", "U1-", "U1-", "U…
#> $ country_code <chr> "Denmark", "Spain", "Spain", "Portugal", "Portugal…
#> $ habitat_code_uri <chr> "http://dd.eionet.europa.eu/vocabulary/art17_2018/…
#> $ habitat_description <chr> "Sandbanks which are slightly covered by sea water…
#> $ geom <MULTISURFACE [m]> MULTISURFACE (POLYGON ((420..., MULTI…Run the following code to have a quick look at the layers geometries
# Transform to Polygon geometry type from Multisurface
if(unique(st_geometry_type(habitats_directive_layers)) == "MULTISURFACE"){
habitats_directive_layers <- habitats_directive_layers %>%
st_cast(to = "GEOMETRYCOLLECTION") %>%
st_collection_extract(type = "POLYGON")
}
# Visualize
map <- mapview(habitats_directive_layers, zcol = "habitat_description", burst = TRUE)
mapFurthermore, we can get data from other EMODnet lots and combine them. The Human Activities portal provides the maritime boundaries of the European Union state members. This time we will not initiate a WFS client, but we will call the service directly. The WFS client will be generated on the fly.
Same as before, we have a look at the layers available first.
emodnet_get_wfs_info(service = "human_activities")
#> ✔ WFS client created succesfully
#> ℹ Service: 'https://ows.emodnet-humanactivities.eu/wfs'
#> ℹ Version: '2.0.0'
#> # A tibble: 95 × 9
#> data_source service_name service_url layer_namespace layer_name title
#> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 emodnet_wfs human_activities https://ows.em… emodnet activelic… Acti…
#> 2 emodnet_wfs human_activities https://ows.em… emodnet aquacultu… Advi…
#> 3 emodnet_wfs human_activities https://ows.em… emodnet baltic Advi…
#> 4 emodnet_wfs human_activities https://ows.em… emodnet blacksea Advi…
#> 5 emodnet_wfs human_activities https://ows.em… emodnet longdista… Advi…
#> 6 emodnet_wfs human_activities https://ows.em… emodnet market Advi…
#> 7 emodnet_wfs human_activities https://ows.em… emodnet mediterra… Advi…
#> 8 emodnet_wfs human_activities https://ows.em… emodnet northsea Advi…
#> 9 emodnet_wfs human_activities https://ows.em… emodnet northwest… Advi…
#> 10 emodnet_wfs human_activities https://ows.em… emodnet outermost… Advi…
#> # … with 85 more rows, and 3 more variables: abstract <chr>, class <chr>,
#> # format <chr>The layer_name for the maritime boundaries seems to be maritimebnds. This dataset was developed based on the official data provided by the European Environmental Agency and the Maritime Boundaries Database compiled by MarineRegions.org (Flanders Marine Institute, 2019).
The sitename variable specifies the type of boundary each feature represents. For illustration purposes, we will filter our request to return only Territorial Seas.
maritime_boundaries <- emodnet_get_layers(service = "human_activities",
layers = "maritimebnds",
reduce_layers = TRUE,
cql_filter = "sitename='Territory sea (12 nm)'")
#> ✔ WFS client created succesfully
#> ℹ Service: 'https://ows.emodnet-humanactivities.eu/wfs'
#> ℹ Version: '2.0.0'
glimpse(maritime_boundaries)
#> Rows: 64
#> Columns: 13
#> $ gml_id <chr> "maritimebnds.49", "maritimebnds.58", "maritimebnds.59", "m…
#> $ objectid <dbl> 49, 58, 59, 60, 61, 62, 54, 55, 56, 57, 39, 40, 41, 42, 44,…
#> $ mblszotpid <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ localid <dbl> 49024, 49064, 49065, 49066, 49087, 49099, 49034, 49036, 490…
#> $ sitename <chr> "Territory sea (12 nm)", "Territory sea (12 nm)", "Territor…
#> $ legalfound <date> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
#> $ legalfou_1 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
#> $ country <chr> "Latvia", "France", "France", "France", "Spain", "Cyprus", …
#> $ nationalle <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
#> $ nutscode <chr> "LV", "FR", "FR", "FR", "ES", "CY", "SE", "UK", "FR", "FR",…
#> $ mblsds_mbl <chr> "In www.marineregions.org", "In www.marineregions.org", "In…
#> $ shape_leng <dbl> 12.005746, 1.949563, 2.072975, 1.837879, 27.878496, 13.6428…
#> $ the_geom <MULTICURVE [°]> MULTICURVE (LINESTRING (23...., MULTICURVE (LINE…Add the maritime boundaries to the previous map with this line. To do so, cats to MULTILINESTRING geometry type, transform to the same crs as the habitats directive sf (only required because of the next step), and crop using the bounding box of the habitats directive sf.